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Abstract 


Numerical investigation of flow and heat transfer in a rectangular duct with 
a built-in circular tube has been carried out for a Reynolds number of 1000 
and blockage ratio of 0.44, A finite volume discretization method due to 
Eswaran and Prakash (1998) and a SMAC (Harlow and Amsden, 1970) 
based solution algorithm have been applied to solve the governing 
differential equations. Since the heat transfer in the duct is dictated by the 
flow structure, the present study is directed toward characterization of the 
flow structure. To this end, the topological theory shows the promise of 
becoming a powerful tool for the study of the flow structure. Computations 
show the helical vortex tubes in the wake and existence of the horseshoe 
vortex system. The w component of velocity is surprisingly very large in 
front and in the near wake of the tube. The qualitative analyses of the 
limiting streamlines on the tube and the bottom-plate reveal the existence of 
a complex flow field. The separation lines as well as singularity points 
(saddle and nodal points) have been investigated. The iso-Nusselt number 
contours and the span-averaged Nusselt number in the flow passage draws 
an insightful conclusion about the heat transfer performance in the duct. 
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Chapter 1 
Introduction 

1.1 Description of the Problem 

Fin-tubes are commonly used in the gas-liquid crossflow heat exchangers, 
where the gas generally flows across the tubes and the liquid flows inside the 
tubes. Figure 1.1 shows a schematic diagram of the core region of a fin-tube 
heat exchanger. The purpose of the fin is to enhance the heat transfer rate on 
the gas side, since the transport coefficient on the gas side is usually smaller 
than the liquid side. Numerical investigation on the related topic has been 
reported earlier by Biswas, Mitra and Fiebig (1994). In the above 
investigation the enhancement of heat transfer from the fin surfaces is 
achieved by inducing longitudinal streamwise vortices in the flow field. The 
longitudinal vortices are generated by placing delta-winglet type vortex 
generators on the flat surface (Figure 1.2). The longitudinal vortices 
developed along the side edge of the delta-winglets due to the pressure 
difference between the front surface facing the flow and the back surface. 
The longitudinal vortices are so called the streamwise vortices since the axes 
of the vortices are aligned to the flow direction. The streamwise vortices 
interact with an otherwise two-dimensional boundary layer and produce a 
three-dimensional swirling flow that mixes near wall fluid with the free 
stream. This mechanism strongly enhances the entrainment of fluid from the 
periphery to the core region of the flow field. Thus the thermal boundary 
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Liquid 
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Schematic Diagram of Core Region of a Fin-Tube Heat Exchanger. 

Figure 1.1 



Winglet pair 


Delta-Winglet Type Vortex Generators on Flat Surface 

Figure 1 .2 




layer is disrupted and the heat transfer rate is enhanced. The additional 
pressure losses are modest because the form drag for such winglet-type 
slender bodies is low. 

The air-cooled condenser of geothermal plants consists basically of the same 
fin-tube arrangement as shown in Figure 1.1. Air is forced through several 
rows of this fin tubes by large fans. The condenser units can be very large, 
consuming a large fraction of the overall capital cost of these plants. In 
addition the power required to operate the fans represents a significant 
parasitic house load, reducing the net power production of the plant. In order 
to analyze the mechanisms involved in the heat transfer and flow behavior in 
such heat exchangers, a detailed investigation on a heat exchanger module is 
necessary. Such a module is shown in Figure 1.3. 

1.2 The Scope of the Present Work 

The main aim of the research is to prescribe the heat transfer enhancement 
strategies of the air-cooled condensers for the geothermal plants. In order to 
achieve this desired objective a detailed three-dimensional numerical model 
has been formulated to provide a better understanding of the flow physics. In 
the numerical model, the full Navier-Stokes equations together with 
governing equations of energy are solved in a rectangular channel with a 
built in circular cylinder (Figure 1,3). A detailed analysis of the flow 
structure along with heat transfer characteristics in such a module is studied, 
as it is crucial for design of high-performance heat exchangers. 
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Heat Exchanger Module 
Figure 1 .3 



1.3 Layout of the Thesis 


In Chapter- 1, we have already discussed the genesis of the problem. 
Chapter-2 of the thesis provides a review of literature relevant to the basic 
mechanisms involved in augmentation of heat transfer and the algorithms 
related to solving the Navier-Stokes equations. A concise review of 
modeling aspects of turbulent flow has also been focussed in this chapter. 
The mathematical formulation of the problem for simulation is presented in 
chapter-3. In this chapter, the geometry, the governing equations, boundary 
conditions and the grid generation techniques are described. Chapter-4 deals 
with the solution algorithm. In this chapter, the discretization schemes and 
the solution procedure using control volume formulation are described in 
detail. Chapter-5 discusses results for laminar flow. In this chapter, a 
detailed study of flow physics and associated heat transfer characteristics are 
presented. Chapter-6 includes the concluding remarks and the scope for 
further research. 
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Chapter 2 


Literature Survey 

2.1 Introduction 

The objectives of present study have been briefly outlined in the previous 
chapter. A variety of experimental and considerable amount of analytical 
and computational research have been carried out on the enhancement of 
heat transfer. Especially the enhancement of heat transfer through 
manipulation of surface geometry has been analyzed by many researchers 
and practitioners since the earliest documented studies of heat transfer. In 
this chapter a survey of the relevant literature is presented to indicate the 
extent of work already reported in open literature pertaining to the 
enhancement of heat transfer by using surface-mounted protrusions. In order 
to analyze the flow structure and heat transfer in such applications, a detailed 
computational study is needed. The flow regime for such applications can be 
both laminar and turbulent. In the present literature review the attention is 
also focussed, on important numerical techniques developed in the recent 
past for solving the conservation equation for complex flows and heat 
transfer. Although the numerical solution of the governing conservation 
equations for laminar flows has become a realizable goal, the solution for 
turbulent flows even for simple geometry is indeed a formidable task. The 
relevant literature on turbulence modeling has been looked into in brief with 
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the intent to understand various modeling and computational strategies. This 
survey helps in suggesting the work that should be carried out to accomplish 
the objectives enumerated in the earlier chapter. The literature is reviewed 
from three different viewpoints. In the first section, an overview of the work 
done by various researchers in the area of augmentation of heat transfer is 
presented and discussed. The second section presents a review of 
investigations pertaining to the numerical methods with regard to the task of 
computing flow fields in complex geometries. Basically, this section deals 
with available scheme for solving the complete Navier-stokes equations. In 
the third section a brief review on modeling aspects of turbulent flow has 
been presented. 


2.2 Fin-Tube Heat Exchangers 


As it was mentioned in the previous chapter that our primary interest is to 
enhance the heat transfer in the gas side of fin-tube heat exchangers (with 
flat fins) and enhancement of heat transfer in fin-plate heat exchangers. With 
this intent, we would like to study different investigations related to 
augmentation of heat transfer suitable for the applications enumerated 
earlier. 


Augmentation of heat transfer is of special interest in channel flows where 
the rate of heat transfer between the fluid and channel walls deteriorates as 
the boimdary layer grows on the channel walls and the flow tends to become 
fully developed. Protrusions can be mounted Qn these channel walls in order 
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to disturb the growth of boundary layer and thereby enhance the heat 
transfer between the flowing fluid and channel walls. Two relevant 
applications using such flow configurations are the heat transfer between the 
gas and the fin in the case of gas-liquid fin-tube crossflow heat exchangers 
and the heat transfer between the flowing fluid and plates in the case of fin- 
plate heat exchangers. The evolution towards a fully developed flow can be 
disturbed by using a multi-louvered surface geometry for plates. 
Investigation by Achaichia and Cowell (1988) provides a detailed 
performance data for louvered fin surfaces. However in using louvered fins, 
enhancement is obtained at the price of high-pressure drop. To circumvent 
this difficulty, protrusions in the form of slender delta-wings or winglets can 
be deployed (Figure 1.2). As shown, the base of the wing remains attached 
to the fin and the apex faces the incoming stream with an angle of attack 
with this configuration. The longitudinal vortices are generated along the 
side edge of the wing-shaped vortex generator due to the pressure difference 
between the fi-ont surface facing the flow and the back surface. These 
longitudinal vortices, generated by the vortex generators, can be made to 
disturb the growth of boimdary layer in a channel by exchanging the fluid 
from the near-wall-region with the channel-core-region and thus they can 
serve to enhance the heat transfer rate while producing less increase in 
pressure penalty. Use of longitudinal vortices for boundary control is well 
known (Pearcey, 1961) and the vortex generators are used in commercial 
airplanes for this purpose. 
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Formation of streamwise longitudinal vortices behind a slender aerodynamic 
object is a research topic of considerable interest for several years. Both the 
theoretical and experimental investigations on flow past a delta wing have 
been conducted and reported in literature by number of researchers. 
Hummel and Srinivasan (1967) have made important contribution in 
revealing the complex flow structure behind a delta wing. They have 
conducted experiments and presented pressure distribution and vortex 
structure of the flow around a delta wing of unit aspect ratio with an angle of 
attack of 20°. 

Thomas et al. (1990) have computed low-speed laminar flow over a low 
aspect ratio delta wing up to 40° angle of attack using an upwind biased 
finite volume algorithm. The differencing schemes used are second-order 
accurate and a multigrid algorithm is employed to promote convergence to 
steady state. The predicted results and lift coefficient have remarkable 
agreement with experiments due to Hummel (1973). At 40° angle of attack, 
a bubble-type vortex breakdown is evident in the computations. 

Experimental investigations due to Fiebig et al. (1986), Fiebig et al. (1991) 
and Tiggelbeck et al. (1992) can be referred to in coimection with 
augmentation of heat transfer by means of longitudinal vortices. 
Experimental investigation due to Fiebig et al. (1986) is the first systematic 
study to compare the performance of different kinds of vortex generators 
viz., delta wing, rectangular wing, delta winglet pair and rectangular winglet 
pair in the Reynolds number range of 1360 and 2270. Their observation 
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depicts that the delta wing is the best vortex generator from the heat transfer 
point of view. Another important feature of their observation is that the heat 
transfer coefficient increases with the increase in angle of attack till the 
vortex breakdown takes place, Tiggelbeck et al. (1992) used multiple rows 
of vortex generators in an aligned arrangement within a channel and 
observed their influence on flow structure. The flow structure in the wake of 
the second row is qualitatively similar to that of the first row. Their flow 
visualization by a laser light shift technique has revealed that the 
concentrated vortex pair generated by a small aspect ratio delta wing at large 
angle of attack has elliptic shape due to the influence of channel walls. They 
also observed that the peak value of the spanwise-averaged Nusselt number 
at the wake of second row is strongly dependent on the spacing between the 
two rows. 

Computational studies on related topics have been performed by Fiebig et al, 
(1989) and Biswas and Chattopadhyay (1992) for laminar flow in 
geometrical configuration of delta wing placed inside a channel. Both 
studies have discussed the influence of angle of attack and Reynolds number 
on velocity and temperature fields. 

It has been already discussed that the gas side heat transfer coefficient in a 
gas-liquid fin-tube crossflow heat exchanger is small compared to that of 
liquid side. The mechanism of heat transfer between the gas and the solid 
surfaces in such cases is to be understood in detail. Let us consider only one 
tube in a channel, formed by the two neighboring fins. The flow field 
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consists of a horseshoe vortex system, a dead water zone at the jimcture of 
the tube and the plate and a von-Karman vortex street in the middle. In order 
to enhance heat transfer in such flow configuration, vortex generators can be 
mounted on the plates, which makes the flow field extremely complex. Dong 
(1989) and Valencia (1992) have conducted experimental investigation, to 
observe the influence of winglet type vortex generators in a channel with a 
built-in circular tube. Biswas et al. (1994) have observed that in the absence 
of the winglet-type vortex generators, relatively little heat transfer takes 
place in the downstream of the circular tube which is the recirculation region 
with low velocity. However, they observed an enhancement of heat transfer 
as high as 240 percent in the wake region behind the cylinder in the presence 
of winglet-type longitudinal vortex generators. Fiebig et al. (1994) have 
studied in detail the conjugate heat transfer of a fin-tube heat exchanger for 
three-dimensional thermally and hydrodynamically developing laminar 
flows. They have identified an unique heat transfer phenomenon, a 
directional reversal of heat transfer (HTR) occurred locally on the fin in the 
tube wake for large Reynolds number and a small fin efficiency parameter, 
Fi (ratio of fm to fluid conductivity times fin thickness to fm pitch). 

Recently computational studies have been extended to periodic 
configurations, which are realistic in terms of plate-fin exchangers with 
perforated fins or slotted channels. Heat transfer surfaces periodically 
interrupted along the streamwise direction have been widely used in order to 
obtain improved performance of the heat exchanger devices. Such an 
arrangement maybe viewed as a succession of collinear plate-segments 
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aligned parallel to flow with slots between successive plates. Each slot 
enables the interruption of thermal boundary layer where the highest 
resistance to heat flux occurs and a new boundary layer with lower thermal 
resistance begins. It also allows a very good mixing of flow with self- 
sustained oscillation when the system is operating above the critical 
conditions (Reynolds number). The channels, which have such operating 
conditions, have been classified as communicating channels. As it has been 
said, above a critical value of Reynolds number (much within the laminar 
region), the initial steady flow becomes unstable to small disturbances. At 
this instant the flow is said to have undergone a bifurcation from steady state 
to time periodic, and finally to self-sustaining oscillatory regime. 
Supercritical bifurcation and narrow band excitation characterize these flows 
(Ghaddar et al.l986 a; Ghaddar et al.l986 b; Amon and Patera 1989). 
Appreciable transport augmentation takes place in this range of Reynolds 
number (Majumdar and Amon, 1992). 

Webb and Ramadhyani (1985) have analyzed fluid flow and heat transfer 
characteristics through a parallel plate channel with transverse ribs, which 
act as vortex generators. Significant heat transfer augmentation is obtained 
in their study. Acharya et al. (1991) explained vortex interaction between 
two such ribs in a channel for the sub harmonic excitation of shear layers. In 
a similar way a significant augmentation is obtained by positioning vortex 
generators at key locations in the flow (Myrum et al. 1992). 
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Many researchers have observed longitudinal vortices in various complex 
flow configurations. Naturally, the interaction of such vortices with 
boimdary layer and its effect on heat transfer is a subject of interest to them. 
Some examples are Taylor-Gortler vortices in boundary layers on concave 
curved surfaces, the horseshoe vortices formed by an obstruction protruding 
from surface and wingtip vortices impinging on a downstream surface. The 
embedded vortex is capable of strongly perturbing the boundary layer and 
influencing the heat transfer characteristics. In addition, longitudinal vortices 
usually maintain their coherence over a long streamwise distance. As a 
consequence, the heat transfer effects behind the vortex generator are very 
persistent. Eibeck and Eaton (1987), Westphal and Mehta (1987) have 
worked extensively in this field with more focus to the turbulent boundary 
layer. Eibeck and Eaton (1987) have conducted experiments on longitudinal 
vortices embedded in a turbulent boundary layer and resultant heat transfer 
effects. Longitudinal vortices are found to influence heat transfer behavior 
significantly. Local Stanton number is increased is by as much as 24 percent 
resulting in a net increase in spanwise average heat transfer coefficient. 
Despite the presence of turbulent diffusion, the influence of longitudinal 
vortices on momentum and energy transport can be traced at a location as far 
downstream as 60 wing chords behind the delta wing. 
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2.3 Numerical Methods for Solving Navier- 
Stokes Equations 


Different solution techniques for full Navier-Stokes equations have been 
developed during the past three decades. The major difficulty encountered 
during the solution of incompressible flows arises from the absence of any 
explicit equation for pressure and due to the nature of spatial coupling of the 
pressure and velocity. For incompressible flow problems, pressure does not 
have the usual thermodynamical meaning. Here it is a relative variable, 
which adjusts itself instantaneously for the condition of zero mass 
divergence to be satisfied at all computational cells. This behavior is related 
to the well-known fact that in an incompressible fluid the speed of the sound 
is infinite. As a consequence, the pressure field cannot be calculated by an 
explicit time advancement procedure, instead it requires at least a partially 
implicit determination which takes into account the coupling between the 
pressure and velocity fields as well as the effects the velocity boundary 
conditions. This aspect is the most distinctive feature of the primitive 
variable formulation of the incompressible Navier-Stokes equations. 

The difficulty of determination of pressure field can be resolved in the 
stream fuction-vorticity approach. But the stream fiiction-vorticity approach 
loses its attractiveness when the three-dimensional flow is computed because 
of the absence of a single scalar stream function in the three-dimensional 
space. 
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A three-dimensional problem demands a primitive variable approach. Efforts 
have been made so that two-dimensional as well as three-dimensional 
problems could be computed following a primitive variable approach 
without encountering non-physical wiggles in pressure distribution. As a 
remedy, it has been suggested to employ a different grid for each variable. 
Harlow and Welch (1965) have used such a staggered grid for the dependent 
variables in their well-known MAC (Marker and Cell) method. The MAC 
method of Harlow and Welch is one of the earliest and widely used explicit 
methods for solving the full Navier-Stokes equations. In this method, 
solutions of velocities are obtained using a two step procedure. In the first 
step called the predictor step, the provisional values of velocity components 
are computed explicitly using advection, diffusion and pressure gradients of 
the earlier time steps. This explicitly advanced provisional velocity field 
may not necessarily ensure a divergence free velocity field. Hence, in the 
second step called the corrector step, the pressure and velocity components 
are corrected so that the velocity field satisfies the continuity equation. This 
is done through the solution of a Poisson equation for pressure. A relative 
technique, known as pseudo-compressibility method, developed by Chorin 
(1967) involves a simultaneous iteration of the pressure and velocity 
components. Vicelli (1971) has shown that the method due to Chorin (1967) 
and the MAC method are equivalent. 

The original version of MAC method has been modified by Harlow and 
Amsden, popularly known as Simplified MAC (SMAC, 1970), Nicholas and 
Hirt (1971) and Hirt and Cook (1972) for application of free surface flows. 
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The MAC method uses a layer of imaginary cells around the boundary of the 
physical domain necessitating the updating of the boundary conditions after 
every change in internal velocity and pressure values. The MAC method has 
been extensively used by many researchers to solve for the flow field in 
complex geometries. For example, Braza, Chassaing and Ha-Minh (1986)' 
have obtained the unsteady wake behind a circular cylinder using the MAC 
method. Mukhopadhyay, Biswas and Sundararajan (1993) have obtained the 
periodic wake behind a rectangular obstacle. In fact, MAC method has been 
successfully used even to simulate highly unsteady turbulent flows 
(Robichaux, Tafti and Vanka, 1992). It has been experienced that the MAC 
method is indeed very efficient in study of temporal flow development but it 
has stability restrictions on the time increment, which slows down the 
calculation for steady flow considerably. 

Since implicit methods have no such restrictions, they are more attractive. 
Patankar and Spalding (1972) have introduced an efficient method know as 
SIMPLE (Semi-Implicit Method for Pressure Linked Equations). This 
method is based on a finite-volume discretization of the governing equations 
on a staggered grid. In order to improve the convergence involved in the 
pressure-velocity coupling, several variants of SIMPLE algorithm have been 
developed. The SIMPLER algorithm of Patankar (1981) and the SIMPLEC 
algorithm of Van Doormaal and Raithby (1984) are improvements on 
SIMPLE. Although the changes to incorporate SIMPLEC into SIMPLE 
algorithm are minor, the consequences can be great as it eliminates the 
approximations made in SIMPLE while deriving the pressure-velocity 
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corrections. A comparative illustration of the operator splitting algorithm, 
viz., the PISO of Issa (1986) and the SIMPLE family of algorithms have 
been reported by Jang, Jetli and Acharya (1986). 

Garg and Maji (1987) have applied SIMPLEC (Van Doormaal and Raithby, 
1984) method for solution of viscous flows through periodically converging 
-diverging tubes. In another numerical investigation, Velusamy and Garg 
(1993) have solved the complete set of Navier-Stokes equations for three- 
dimensional developing flow in elliptic cross-section ducts. Application of 
finite volume methods using non-orthogonal coordinates and collocated 
grids is reported by Rhie and Chow (1983) and Peric (1985). Peric et al. 
(1988) have shown that collocated arrangement converges faster than the 
staggered variable arrangement and has advantages when extensions such as 
multigrid techniques and non-orthogonal grids are considered. It maybe 
worthwhile to note from Benard and Thompson (1984) that while staggered 
grid is a direct remedy to avoid pressure split with conventional finite 
difference interpolation of flow variables, proper interpolation of velocities 
still remain crucial to avoid oscillatory velocity field. An appropriate 
pressure interpolation avoids the pressure splits in collocated grid 
arrangement. A finite volume based procedure has been successfully 
developed using such collocated velocities and pressure. The concept of 
momentum interpolation of cell face pressure from nodal values has 
established an effective equivalence in the methods for staggered and 
collocated arrangements. Mukhopadhyay et al. (1993) has developed a 
numerical method for predicting viscous flows for incompressible 
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geometries. Integral mass and momentum conservation equations are 
deployed and these are discretized into algebraic form through numerical 
quadrature. The physical domain is divided into a number of non-orthogonal 
control volumes which are isoparametrically mapped on to standard 
rectangular cells. Numerical integration for unsteady momentum equation is 
performed over such non-orthogonal cells. The algorithm is tested on some 
complex problems. The results exhibit good accuracy and justify the 
applicability of the algorithm. 

Kobayashi and Pereira (1991) have modified the momentum interpolation 
method suggested by Peric (1995) and named it as Pressure-Weighted 
Interpolation Method-Corrected (PWIMC). In this method, the non- 
orthogonal terms in the momentum equations were solved explicitly, 
whereas in the pressure-corrections they were dropped. 

It is clear that substantial progress has been made for the development of 
algorithms for complex incompressible flow simulations but none of the 
prescriptions is universal. Depending on the nature of flow and geometry 
etc., one can always go for the best-suited discretization procedure (see 
Vanka, 1987; Fletcher, 1988). It is very difficult to be conclusive about any 
algorithm with its universal sense of applicability. However the finite 
volume algorithm developed by Eswaran and Prakash (1998) for solving 
time dependent three-dimensional full Navier-Stokes equation has the 
following advantages: 


18 



• It uses collocated grid arrangement, hence all the advantages of 
collocated grid arrangement are gainfully utilized. 

• The essence of the algorithm is SMAC so it is suitable for modeling 
unsteady flows (see Kim and Benson, 1992). 

• The governing equations are discretized in the physical plane itself 
without coordinate transformation. 

This algorithm has been used to solve the flow problem and extended to 

solve the energy equation for the present work. 



Chapter 3 


Mathematical Formulation 

3.1 Introduction 

Arising out of the practical importance of the flow configuration and thermal 
conditions as described in the previous chapter, we a need is felt to compute 
three-dimensional flow in a horizontal channel with a built-in circular tube. 
This numerical simulation is capable of providing with quantitative 
information of flow structure and temperature distribution in the entire 
domain with their dependence on various governing parameters. 

3.2 Statement of the Problem 

The computational domain is shown in Figure 3.1. Two flat plate fins form a 
channel of height H, width B = 11.25H and length L = 25H. The circular 
tube of diameter D = 5H is located at a distance Li= 1.25D from the inlet. 
The blockage ratio, D/B is kept fixed to 0.44. 

The flow field around the circular tube in a channel is characterized by the 
following parameters (Figure 3.1): 

• The Reynolds number. 
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The grid-system and the computational domain 

Fioure 3.1 



• The inverse blockage ratio D/B 

• The position of the cylinder in the channel (Ll/L) 

• The velocity profile at the channel inlet 

• The relative channel height (H/B) 

These parameters determine the structure of the vortical motion in the wake. 


3.3 Governing Equations 

The three-dimensional Navier-Stokes equations for laminar flow of an 
arbitrary spatial control volume V bounded by a closed surface S can be 
expressed in the following general convection-diffusion-source integral 
form: 


d_ 

dt 



V 




J p u • dS = 0 
s 


( 3 . 1 ) 


+ |[ou«i-r^v4dS =Js^dV 

vs V 


Where p represents the fluid density, u is the fluid velocity, (f) stands for 
any vector component or scalar quantity, is the volumetric source term. 
For incompressible flow of Newtonian fluid, the equation takes the form 
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(3.3) 


Ju.dS=0 

s 


_a 

dt 


^p<j>iY + J[pu«> - V«i].dS = fs^dV (3.4) 

vs V 

and the source term for momentum equation becomes — fp I.dS where 

P J 


is the unit tensor. 


In this formulation we work with Cartesian components of velocity. So ^ 

can be the three-cartesian components of velocity u, v, w as well as any 
scalar e.g., temperature which need to be determined. It is to be noted that 
equation 3.2 is a perfectly general equation from which the three 
fundamental equations of fluid mechanics (i.e., Continuity, Momentum and 
Energy) can be derived. 

Table 3.1 : Values of the variables in the general transport equation 3.2. 


Equation 




Continuity 

1 

0 

0 

Momentum 



ap 

Sxj 

Energy 

T 

k/cp 

0 
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3.4 The Grid 


Figure 3.1 shows the three-dimensional grid used for our present 
computation. The grid generation techniques and its properties are described 
in brief in the following sections. 


3.4.1 Grid Generation Technique 

A differential equation method has been used for the generation of the three- 
dimensional grid. The initial grid is generated by the method of transfinite 
interpolation. This method essentially uses linear interpolation scheme to 
compute the interior points by using the coordinate values from the 
boundaries. The algebraic grid so generated in general is not smooth. This is 
because of the fact that the slope discontinuities at the boundary propagate in 
the interior of the domain. So the mesh obtained by algebraic mapping is 
further improved by the use of Partial Differential Equations (PDE). In this 
technique, a system of PDEs is solved for obtaining the location of the grid 
points in the physical space, whereas in the computational space is of 
rectangular shape with uniform spacing. 

All the three types of partial differential equations, namely elliptic, parabolic 
and hyperbolic are used for grid generation. Elliptic equations are common 
for closed geometries. Hyperbolic equations are used for domains where 
outer boundary is not prescribed. Use of parabolic equations provides the 
advantage of both, speed (marching nature, hyperbolic) and smoothness 



(diffusive nature, elliptic). Parabolic grid generators are in development 
stage and many issues like perfection, robustness of technique, specification 
of the control function, are yet to be resolved. The present grid is developed 
on the basis of elliptic grid generators. 


Any grid generator system demands geometric information to be fetched 
from the boundaries. Hence the steady state boundary valued nature of 
elliptic equations makes them most favorite and most widely used grid 
generators. A Laplace equation or a Poisson equation with a Dirichlet 
boundary condition can be used for this purpose. Standard Poisson equation 
is of the form 

V^^=P 

V^77=Q (3.5) 

Where P and Q are known as control functions. 

If P = Q = 0 then 

The equation gets modified to Laplace form as 


v^^ = o 

V^TJ = 0 


(3.6) 


The above equations can be solved by finite difference techniques to get the 
location of the interior points from the boundary. The initial guess needed to 
solve the Laplace equation can be generated using algebraic mapping. 
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3.4.2 A Note on Control Functions 


The laplacian operator can provide highly smooth grids. The grid lines tend 
to be equally spaced in the absence of boundary curvature but they show a 
tendency to concentrate near the concave boundaries. This may prove as a 
drawback when the accumulation of grid lines is enforced in some specific 
region. Such a control on the grid lines can be obtained by proper selection 
ofcontrol functions P and Q. The other objective of the use of control 
functions is enforcing the orthogonality over the selected edges. Such 
boundary orthogonality simplifies the application of boundary conditions. 
The use of control functions is however associated with additional 
complexity in the transformed domain. 

3.4.3 Computational Procedure for Grid 
Generation 

The algorithm used for the generation of the present grid can be outlined as 
follows: 

1 . Take the input in the form of geometric data. 

2. Define the computational grid (on ^ ,77 plane) based on the number 

of grid points. 

3. Define the boundary points along the edges of the physical 
domain. 
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4. Generate an initial guess for the Laplacian by algebraic grid 
generation method. In this step just a linear interpolation scheme is 
used to compute the interior points. 

5. The guess is improved a little by adopting simple four-point 
formula. 

6. The actual Laplacian is in the form of 


ax^^ - 2bx^^ + cx^^ = 0 
- 2by^^ + cy^^ = 0 

where 

2 2 
a = x„ + y„ 

rf J T] 

b = x^x^ + y^y^ 
c = x]+y] 

is solved till convergence is obtained to get the final grid. 


(3.7) 


3.4.4 Grid Property Evaluation 

Truncation error of numerical solutions depends upon the grid on which the 
discretized approximation of the differential equation is solved. A desirable 
grid should have minimum truncation error. Certain properties of the grid 
are found to have strong effect on the truncation error. Parameters that effect 
these properties can be used to express the quality of the grid. The 
parameters are 
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• Transformation Jacobian. 

• Skewness. 

• Aspect ratio. 

• Adjacent cell ratio. 

1. Transformation Jacobian 


The transformation Jacobian represents the local area-scaling factor for the 
grid generated. The determinant of the Jacobian matrix should not become 
zero if the mapping is to remain one to one. The quantity J =0 will lead to 

a situation of an elemental area (or volume) in the transformed domain being 
mapped onto a point in the physical domain. Such a mapping is not 
acceptable. 


2. Skewness 


Skewness is a measure of non-orthogonality between two intersecting grid 
lines. Higher the skewness, higher will be the truncation error due to loss of 
independence between the grid lines. Many methods are available to 
evaluate the cell skewness. These methods are stated in the following section 

• Intersection Angle 

The intersection angle, 0 between two grid lines is calculated. It is desirable 
that the intersection angle be 90°, however, this is not always possible. At 
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least, it is desirable to maintain the intersection angle in the range of 
45°< e< 135°(Figure3.2) 


• Length of diagonals 

In this method the skewness is calculated by using the length for the 
diagonals of a grid cell. For the grid cell ABCD as shown in Figure 3.3, the 
skewness is calculated as 


Skewness = 1 - 


Length of smaller diagonal 
Length of 1 arger diagonal 


BD 

AC 


• Area of triangles formed by the diagonals 

In this method, skewness is calculated using the area of the triangles formed 
by the diagonals and the four sides of a grid cell. 

Area of the smallest triangle AODC 

Skewness = = 

Area of the biggest triangle A O AB 

In the above two methods, normally the maximum allowable limit for 
skewness is taken around 0.5. 

• Aspect Ratio 

This is the maximum ratio between any of the two adjacent sides of a cell. 
Aspect ratio of less than 6 is desirable. It should not exceed 20. 
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• Adjacent Cell Ratio 

This is the measure of grid smoothness. It is the ratio of the areas of two 
adjacent cells. It is desirable to have it close to one. 


3.5 Boundary Conditions 

The governing differential equations are elliptic in space and parabolic in 
time. We need the boundary conditions for all the confining surfaces. The 
boundary conditions of interest in the present investigation are 

• Top and bottom plates 

u = V = w = 0 (No Slip boundary condition) 



dz 


T = T^ (T^ represents wall temperature) 


• Side Wall 


dn _ dw 
~dy~'~dy 


V = 0 (Free Slip boundary condition) 
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ay 


= 0 


dT 

dy 


= 0 


• Channel Inlet 

u = Uoo, v = w = 0 


ax 


= 0 


T = T 

X Xqq 


• Channel Exit 

There is no unique prescription for outflow conditions. The evolution of 
mass flux through open boundaries with time is found by means of a discrete 
radiation condition (Orlanski 1976), which allows changes inside the flow 
field to be transmitted outward, but not vice-versa. Thus, a disturbance (such 
as a moving eddy) created inside the computational domain can leave the 
domain through an open boundary, but it cannot reenter. The discrete 
radiation condition allows the calculations to be made with fewer grid cells 
(i.e., less storage and computer time) than might otherwise be necessary with 
fixed flux condition. Orlanski condition is expressed as follows: 
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^ dcf) 

+ Uav. — = 0 (where ^ represents u, v, w or T) 


P = Poo 

• Obstacle 

u = V = w = 0 


^ = 0 
an 


(where n signifies normal direction) 


T = T 
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Chapter 4 


Discretization Procedure and Solution 
Algorithm 

4.1 Introduction 

In order to solve numerically the governing equations in the prescribed 
computational domain with given initial and boundary conditions, two major 
discretization steps have to be performed. First the domain must be 
discretized, resulting in the numerical grid used for the computation. An 
explanation about the grid generation procedure was given in Chapter-3. The 
governing equations are discretized on the basis of this numerical grid. 
Various ways of discretization can be applied to the above set of partial 
differential equations. The most common methods are those using finite 
differences, finite volumes and finite elements. In the present work, a finite 
volume scheme was chosen for the spatial discretization of the governing 
equations. Thereby, a conservative formulation of the discretization was 
possible, assuring automatically the satisfaction of the global balance of the 
conserved quantities, independently of the coarseness of the numerical grid. 
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4.2 Finite-volume Method 


The conservation equations are discretized employing the finite- volume 
approach of Eswaran and Prakash (1998). The solution domain is divided 
into a number of contiguous (finite) control volumes (CV). The control 
volumes are defined by the coordinates of their vertices, which are assumed 
to be connected by straight lines. The coordinates of the control volume 
vertices are calculated by the grid generation procedure. A collocated grid 
arrangement is employed and all the dependent variables (u, v, w, p, T) are 
defined at the same location, the centroid of the control volume (Figure 4.1). 
E, W, N, S, T and B indicate the six neighboring control volume centers for 
the east, west, north, south, top and bottom neighbors. The face center points 
e, w, n, s, t and b are located at the intersection of the lines joining the 
midpoints of the opposite edges. For example, te, be, ne and se are the 
midpoints of the edges that form the east face with e as the center of the cell 
face. 

The collocated grids have the following advantages over the staggered grids, 
especially when the non-orthogonal coordinates are used: 

• All variables share the same location; hence, there is only one set of 
control volume. 

• The convection contribution to the coefficients in the discretized 
equations is same for all variables. 
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• For complex geometries Cartesian velocity components can be used in 
conjunction with non-orthogonal coordinates, yielding simpler equations 
than when coordinate oriented velocity components are employed. 

• There are fewer constraints on the numerical grid, since there is no need 
to evaluate the grid sensitive curvature terms. 

Before the time stepping can be done it is necessary to calculate geometrical 
parameters. These parameters include the surface vectors for the six faces e, 
w, n, s, t, b for each finite volume and its volume. 


4.3 Surface areas and Volumes 

The finite volume vertices are numbered 1 to 8 in the manner shown in 
Figure 4.1, The outward surface normals and volume can be found in the 
following maimer as suggested by Kordulla and Vinokur (1983) and 
Eswaran et al.(1995) 

Defining r y = - rj where tj and rj are the position vectors of points i 

and j respectively, we have for 


• East Face 


Se = ^(>'74 X •'Ss) 


(4.1) 
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West Face 


Sw =■^(*'16 X J - Sl ) 


(4.2) 


North Face 


Sn =|(>-27 X «-63) 


(4.3) 


• South Face 


Ss =^(ri8 X r45) 


(4.4) 


Top Face 


St =^(•‘75 X rgg) 


(4.5) 


Bottom Face 


Sb =|('‘13 ^ •■24) 


(4.6) 


The volume of the cell is calculated from the cell coordinates, with the 
assumption that the linear segments to form the six cell faces join the cell 
comers (Eswaran, 1995). 
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(4,7) 


r7i.(S3+Sb+S„) 

4.4 Discretization Procedure 

The main steps of the discretization procedure to calculate the convection 
and the diffusion fluxes and source terms are outlined below. The rates of 
change and source terms are integrated over the cell volume, whereas the 
convection and the diffusion terms form the sum of fluxes through the faces 
of the control volume. 


4.4.1 Discretization of the Continuity Equation 

Equation 3.1 is discretized in the following way: 

j* yOU-dS w ^p(u*S)j = ^puj -Sj (4.8) 

j=e,w,n,s,t,b j 

•th. 

where S j is the surface vector representing the area of the j cell face and 
Uj is the velocity defined at the face center j. 

The discretized form the continuity equation is as the following: 

^Fj=Fe+F^+F„+Fs+Ft + Fb=0 (4.9) 

j 
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where the Fj is the outward mass-flux through face j, defined by; 

Fj = pUj.Sj 

4.4.2. Discretization of the General Equation 

^£^) Rate of Change 

The value of the dependent variable (j) at the centroid of the control volume 
(the geometric center) represents an average over the CV as a whole. Thus 

1 f ( 4 . 10 ) 

dt Jv At At 

where V is the volume of the cell. 

(T)) Convection Fluxes 

The surface integral over convection flux 
in the following form 

j 

where is the value of ^ at the center of face j. Thus 


of variable (j) can be approximated 


( 4 . 11 ) 

j 
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( 4 . 12 ) 


f -dSK 5 + Fn(i„ + Fs(!is + F,(it +Fb«Sb 

JS 

where, is the (interpolated) value of the variable^ at the east face center, 
etc. This can be evaluated by using a central difference linear interpolation 
etween the neighboring nodal values and^Z^ . At east face the value of (f)^ 


is given by 



Vi 


E 


Ve+Vp 


(fy^ + 


Vt 


Ve+Vp 




( 4 . 13 ) 


where Ve and Vp are volumes of the cells aroxmd the points E and P 
respectively and (j> e and^ pare the values of the dependent variables at 
these points. In a collocated grid system, all dependent variables are defined 
at the same location hence exactly the same interpolation scheme is used to 
express all of them at the interfaces. The central difference approximation to 
compute the convection may lead to numerical stability problems. Therefore 
the convection flux is blended with a first-order upwind differencing scheme 
(UDS), and the difference between the central difference scheme (CDS) and 
UDS approximations as 

Fe?Se = (Fe?»e)™® + -(Fe^)™®] (4-14) 

The upwind differencing scheme is based on the assumption that the 
convected cell face value is equal to that at the upstream cell along the same 
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coordinate direction. Thus, the value (j)^ at the east face is assigned the value 
^ if Ug > 0 , i.e., the flux Fg is negative. This can be conveniently 
summarized as 




P [[Fe E Fg 


r 


Ve . , Vp 


-^p+ 




Vg + Vp Vg + Vp 

(4.15) 


The symbol [[,]] signifies greater of the two quantities enclosed inside the 
brackets. Similar expressions can be written for the other cell faces (Eswaran 
etal. 1995) 


Fw^w = ^pttFw»0]]-^^[[-F.^y,0]] + 


• w 


Vr 




V 




w 


l,V^+Vp'^’ Vw+Vp 


^p 


“ ^ F^ ,0]] \ 


(4.16) 


Fni^n = «>pEF„,Oll-«iN[I-Fn,0]] + r 


■n 


Vi 


N 


-^p 


Vr 


V 


Vn+Vp Vn+Vp 




-<^pttFn.0]]+<^Ntt-Fn,0]] 


(4.17) 
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Fs4 = + r iFs 


-^p[lF„0]]+^s[[-Fs,0]] 


r 


Vp 


-^g- 


Vs 


(f>Y 


^Vg+Vp^" Vs+Vp 

(4.18) 


F,A = «>pllFt,0]]-«>TlI-F„Ol]+r 

-«>pttF„0]]+«»TlI-F„Ol] 


( 


Vt ^ , Vp 




-^^P + 




^Vx+Vp^" Vj+Vp 

(4.19) 


Fb<^b = (^p|lFb»^l~<^Btt~Fb>0]]+r ^ 


>p[Fb,0]]+«»B[-Fb.0| 


r 


Vp_x_. Vb 


Vb + Vp 


^B+ 


Vb+Vp 

(4.20) 


^p 




In a fully implicit method the upwind parts of the above equations are 
implicit and they are incorporated in the coefficients of the unknown 
velocity during the pressure-velocity iterations. The CDS terms on the other 
hand are evaluated using the previous iteration values and used as a source 
term on the right side of the same equation. This is the so-called deferred 
correction approach of Khosla and Rubin (1974). Multiplication of the 
explicit part by a factor ofx(0<;'<l) allows the introduction of numerical 

diffusion. For the first order upwind difference scheme, y = 0 and for the 
second order central difference, y = 1. The deferred correction approach 
enhances the diagonal dominance of the coefficient matrix. However, the 
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present solution scheme being explicit, the formal accuracy depends on the 
value of Y . 

(c) Diffusion Fluxes 

Historically Peric (1985) developed the expression for the diffusion fluxes 
for the finite volume approach for the first time. However, the exact 
expressions in three dimensions are not explicitly available in open 
literature. Eswaran and Prakash (1998) have explained the philosophy of 
defining the diffusion fluxes at the cell interfaces in a lucid manner. Eswaran 
and Prakash (1998) evaluate the diffusion flux of a variable (f) through the 
cell faces in the following manner 

-ds » -s). = pf (4.21) 

j=e,w,n,s,t,b j 

For any face we can write 

S j = aiv} + a 2 ^ + cc-i, v? (4.22) 

12 3 

where n , n and n are any three linearly independent (not necessarily 
orthogonal) unit vectors. Therefore 

-Sj = V(l) • ) (4.23) 
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— (j) ' vl' + (f) • 10^ + (f) • 


12 3 

If 5 A^ are the differences in (f) between the two ends of the 

line segments Ax\ Ax^, Ax^,then 

A<2^^ = V<zi-Ax\ A(l>^ =V(f>‘Ax^, A^^=V^-Ax^ (4.24) 


If Ax\ Ax^, Ax^ are in the directions of n\ and respectively, 
then it follows from equation (4.24) that 


A(p^ 

Ax^ 


V^-n\ 


A(f>^ 

Ax^ 




Acfr" 

Ax^ 




(4.25) 


where Ax^, Ax^ and Ax^ are the magnitudes of Ax\ Ax^and Ax^. 
Combining equations (4.23) and (4.25) we have 


V^-S = 


ai 


Ax' 


+ a2 


Ax} 


+ «3 


A(p^ 

Ax^ 


(4.26) 


To obtain ctj, a 2 and a^, we express 
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(4.27) 


n'=(nii ni2 11,3) 
n^=(n2i n22 1133 ) 

n^=(i>3i n32 n33) 


where riu, ni 2 i arc the Cartesian components of n and which can be 

easily determined by where AxJ, Ax 2 , AX 3 are the 

Ax^ Ax"^ Ax^ 

three components of vector Ax^etc. The other values n 2 i...n 33 can be 
similarly obtained. Hence equation (4.22) can be written as 


"nil 

ni2 


T 

“ai“ 


Sij 


1122 

1123 


^2 


S2j 

-^31 

n32 

ll33_ 


_«3_ 


.S3j. 


(4.28) 


where S^j, S 2 j, S 3 j are the Cartesian components of the surface vector 


Using Cramers rule 
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D ’ 




ai = 


«2 = 


0^3 = 


(4.29) 



where D is the determinant of the coefficient matrix. Dj is obtained by 
replacing the first column of D by the column with the elements j , $ 2 j , 


and S 3 j . Thus a^, 0 ^ 2 , and are determined. 

The diffusion flux is made of two distinct parts: normal derivative diffusion 
flux and cross derivative diffusion flux. The second part arises from the non- 
orthogonality of the grid. The normal derivative diffusion flux of ^ through 
any cell face involves the values of (f) at cell centers whereas the cross- 
derivative diffusion flux takes into account the edge center values of (j ) . The 

normal derivative diffusion flux is treated implicitly and is coupled with the 
implicit part of the convective flux to calculate the main coefficients of the 
discretized equations while the cross-derivative diffusion flux is treated 
explicitly to avoid the possibility of producing negative coefficients in an 
implicit treatment. This term together with explicit part of convective flux is 
added to the source term. The example of the east face is taken to illustrate 
the diffusion model (Figure-4.2). Given the edge center values 

^te’ ^se’ <^ne normal diffusion term — , and the 

Ax^ 


cross diffusion term 




e J ^se 
and 




ne 


Ax^ 


Ax^ 


Finally, the diffusion flux is 


computed by 


Fj 


r 




«i — — r- + ^^2 


V 


Ax^ 


Ax^ 


+ 




se 


Ax' 


(4.30) 


/ 
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E 


Face representation to illustrate the diffusion model 

Figure 4.2 



In order to calculate the edge center values for the calculation of cross- 
derivative diffusion flux, the following interpolation scheme is proposed: 

, , Vp V-r , Vp , 

^^te - + 7^: <^TE + 77 ^^E + (4.31) 

''tot ^tot ''tot ''tot 

^tot = '^TE + Vp +V'p -1- Vp 

^he = TT^^^BE + (4-32) 

^tot ''tot vtot ^tot 

Vtot = Vgp -I- Vp -1- Vg + Vp 

VnE Vp Vxt Vp 

Ke = T 7 — <^P + TT^^^NE + TT^^^E + (4-33) 

''tot ''tot ''tot ''tot 

Vtot = Vne + Vp + Vn -hVp 

Vcp Vp Vc Vp 

= 77 — ^P + TT^^^SE + TT^^E + (4-34) 

''tot ''tot ''tot ^tot 

Vtot=VsE+Vp+Vs+Vp 


where, Vpp is the volume of the cell that is located at the top-east of the cell 
P. The . other edge center values of the dependent variables can be 
interpolated in a similar way. 
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(d) Sources 

The source term is to be integrated over the cell volume. The contribution of 
a source term in the discretized equation may be thought of 

J^S^dV»(sJpV (4.35) 

Apart from the real source , explicitly treated parts of the convection and 

diffusion fluxes maybe added to S^. The momentum equations contain 

pressure terms. These terms are also treated explicitly. Its discretization is 

analogous to that of the ordinary diffusion flux, i.e., for the i^^ momentum 
equation the pressure term is 

-JpupS » ^PjSij (4.36) 

J 


fVi fh 

where pj is the pressure at the j face center and Sjj is the i component 
of the surface vector for face j. 


4.5 Pressure-Velocity Coupling 

To obtain the velocity and pressure field satisfying the mass and momentum 
conservation laws, a procedure akin to Simplified Marker and Cell (SMAC) 
method is used. This method offers an efficient and an easy way of pressure- 
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velocity coupling. It is basically a semi-explicit method. The momentum 
equations are discretized in an explicit manner with the exception of the 
pressure gradient terms that are treated implicitly, and the continuity 
equations are also enforced implicitly. This can be expressed mathematically 
by the following two discretized equations of momentum and continuity; 

j j 


and 


= 0 (4.38) 

j 

The momentum equations are solved using the guess values of velocity and 

____ ♦ 

pressure field. The provisional velocity components Uj are calculated from 

the following equation 



Up = 



(4.39) 


where Up is the value of velocity at earlier iteration and is the pressure 

term. This provisional velocity in general will not satisfy the continuity 
equation. The continuity equation (4.38) in another form reads as follows; 


51 


t ■ ■ ■ v.nr? 

6 ? 1 


Ir 




j J 


(4.40) 


where Fj is the uncorrected mass flux obtained from the provisional 

velocities and Fj is the mass flux correction. To evaluate the terms on the 

right side of equation (4.40), it is necessary that the variables (velocities and 
pressure) are to be known at the cell faces. Due to the non-staggered 
arrangement, if the variables at the cell faces are evaluated by linear 
interpolation between the adjacent cell center quantities then the pressure 
iteration does not converge. The solution leads to a Checker Board pressure 
field (i.e., spurious oscillation of pressure may occur), as shown by Rhie 
(1981). This problem can be circumvented by using the concept of 
Momentum Interpolation (Majumdar, 1988). The essence of the concept is 
that the velocities at the cell faces are computed by linear interpolation of the 
convective and diffusive terms but not the pressure term. Thus the method 
with the collocated variable arrangement relies indirectly on the staggering 
idea. Following this idea the interpolated velocity at the east face of the 
control volume given by is obtained in the following manner 
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(4.41) 


Ue = (vp,Ve) Vp 

P 


where 

''p = 

Ve = “e-^(Fe+Fe) (4-43) 

and the overbar indicates a linear interpolation using equation (4.13). So the 
uncorrected mass flux for the east face using equation (4.41) becomes 

F* = pUe-Se = p(vp,VE)-Se - AtVp-Se (4.44) 

Equation (4.44) in its generalized form for any cell face is in the form 

F]" = p V j - S j - AtVp j • S j (4.45) 

To enforce mass conservation, the velocity and pressure corrections are 
introduced, linked by 

u'i Vp'i (4.46) 

P 
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This step simulates the effect of a staggered grid (face-center velocity and 
cell-center pressure) and reaps the benefits of faster convergence. From 
equation (4.46) we get the mass flux correction at the face j 

Fj = pu'j-Sj =-AtVp'j.S'- (4.47) 

Since equation (4.47) is exactly of the same form as equation for diffusion 
flux equation (4.21), with only the variables interchanged, therefore the 
same method can be used for computing Fj 

Combining equations (4.40), (4.45) and (4.47) yields 
-AtVpj-Sj = /?Vj-Sj - AtVpj-Sj 

or AtVp'j-Sj-hSjn= 0 (4.48) 

(S„,=: pVj.Sj -AtVpj-Sj) 


which gives the pressure correction equation. After the solution of the 
pressure correction equation, the nodal velocity, mass fluxes and pressure 

n 4-1 

are updated. The corrected pressure p is obtained by 


P 


n+1 



(4.49) 
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Instead of solving equation (4.37) again with the corrected pressure field, 
velocity calculations can be corrected in the form 



(4.50) 


,1 

The corrected velocity field Uj is then obtained by 

rH" 1 ^ / 

Ui = Uj + u- (4.51) 

Using this new velocity distribution, the energy equation is discretized and 
solved. On integrating the conservative form of the steady state energy 
equation over a control volume (neglecting the dissipative terms as, it is 
significant for high speed flows), we arrive at the following temperature 
equation 


a^(VT.S)j =^(pUj.Sj)Tj 

j j 


(4.52) 


which is solved iteratively to get the temperature distribution of the current 
time step. 

The first terms of equation (4.48) and equation (4.52) are expanded 
following the philosophy that was followed during the discretization of 
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equation (4.21). Subsequent regrouping of the normal derivative terms yields 
the standard Poisson Equation. The Poisson equation can be solved by a 
multitude of methods. We have used the Gauss Seidel iteration method. This 
method is quite reliable and robust. The flow chart with the main steps of the 
algorithm is given in Figure 4.3. 


4.6 Solution Algorithm 


The velocity pressure and temperature fields are calculated by the following 
the procedure of Eswaran and Prakash (1998). The procedure can be 
summarized in the following way 

1. Make initial guess of the pressure velocity and temperature fields. Use 
equations (4.41) and (4.42) to calculate cell center Up and Vp. Use 
linear interpolation (equation 4.12) to obtain the face center quantities 

2. Compute the mass flux through cell face j using equation (4.45) i.e., 

F* = pvj-Sj-AtVpj.Sj 

3. Use equation (4.46) to compute the flux correction at the cell face j using 
equation (4.47) i.e., 
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Figure 4.3 Iterative Solution Scheme 








Fj =-AtVp'j.S^ 


This is computed using the fomulation for diffusion fluxes [equations (4.21 

t 

— 4.34)]. For calculating by Fj , the variable (j) is to be replaced by . 

4. Compute the residue, 5? for each cell, 

j j 

5. Calculate the cell-center pressure correction (using Gauss-Seidel method) 
from the relation 




O) 

Up At 


where O) is the relaxation factor and ap stands for the diagonal coefficients 
of equation (4.48), which is calculated from 


ai 


1 ^2 

«2 


«3 


1 

Ax 

x: 

<1 

1 

e Ax^ 

n Ax^ 

s Ax^ 

s Ax^ 

t_ 


where cti, a 2 etc. are the same as that are used in equations (4.21 - 4.34). 
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6. If 9? rms > ^ S® to step 3. 


7. Store the updated mass flux through cell faces using equation (4.45) 


F* = yOVj -Sj - AtVpj -Sj 


8. Store the updated pressure at the center of the cells 


Pp = Pp + Pp 


9. Store the cell-center corrected velocity 



At 

yoVp 


Zp'iSij 


Up = Up + Up 


10. Calculate the residual for each cell 


s = ;^(puj.Sj)Tj 
j 


a^(VT.S)j 

j 
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1 1 . Calculate the cell-center temperature using Gauss-Seidel method, to 
obtain 


Tp = Tp + 

bp a 

where iQ is the relaxation factor and bp stands for the diagonal coefficients 
of equation (4.52) i.e., calculated from the relation 

bp = — ap 

12. If 5R rms^ ^ 


13. Go to step 1 and repeat the process until the steady state is reached. 
It can be shown that solving equation (4.40) is same as solving 




P- 


V-U; 


At 


and finding the corrections 


t 


u 


At 

P 


Vp' 
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such that 


U; 


n+1 


= U; + U; 


P 


n+l 



+ p 


I 


Thus the procedure above can be called solution of Pressure Correction 
Linked Equations. 


4.7 Numerical Stability Considerations 


Since we are using a Semi-Explicit method, which relies on explicit 
differences, suffer from time step restrictions. For a given mesh the choice 
of the time step is determined through stability analysis which has to take 
care two conditions. First, fluid should not be allowed to cross more than 
one cell in any one time step. This restriction is derived from the Courant- 
Friedrichs-Lewy (CFL) condition given by 


5t 


1 


< mins 


<5x Sy 5z 



(4.52) 
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where the minimum is with respect to every cell in the domain. Typically 
(5^t is chosen equal to one-third to two-third of the minimum cell transient 
time. 


Secondly, when a nonzero value of kinematic viscosity is used momentum 
must not diffuse more than approximately one cell in one time step. A linear 
stability analysis shows that the restrictions on grid Fourier number will 
yield 


1 {SxfjSyfjSz) 

2 (.Jxf + {Syf + {Szf 


( 4 . 53 ) 


Finally, the minimum of the above two time increments is chosen for 
computation. 
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Chapter 5 

Results and Discussion 


5.1 Introduction 

A three-dimensional numerical simulation of the flow structure and the heat 
transfer characteristics in a rectangular channel with a built-in circular tube 
(Figure 1.3) are presented The geometrical configuration represents an 
element of a gas-liquid fin-tube crossflow heat exchanger. The following 
characteristic features of flow are established: 

• The flow that curls around the tube separates at the downstream of the 
leading edge of the tube and forms the recirculation region in the tube 
wake with screw like motion of the helical vortex tube. 

• The wake region of the tube shows strong three-dimensionality. 

• The limiting streamlines on the built-in tube depict typical separation 
lines. 

• The limiting streamlines on the bottom wall reveal presence of th 
saddle point of separation in front of the tube and the bottom wall 
junction. 

• The horseshoe vortices, which form in the front of the tube near the 
bottom and top plates 


The corresponding characteristics of heat transfer influenced by the flow 
structure are also established. 
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• The high heat transfer near the leading edge of the fins and in the 
region influenced by the horseshoe vortices. 

• The poor heat transfer in the wake region of the tube. 

For the computation of the present problem, the finite-volume method of 
Eswaran and Prakash (1998) is used. A 100 x 51 x 17 grid-mesh is used in 

the present computation. The divergence-free criterion is satisfied using an 

upper bound of 10 . The blockage ratio (D/B) and Reynolds number are 

kept fixed to 0.44 and 1000 respectively. In this work no symmetry is 
assumed in the midplane of the channel. 

5.2 Flow Characteristics 

In this section a detailed investigation of the flow pattern in the passage of 
the fin tube element (Figure 1.3) is presented. For a three-dimensional flow 
simulation, it becomes quite difficult to handle the enormous data and only 
effective use of some approaches proposed in literature can depict the flow 
stmcture. Of the possible methodologies that allow us to explore the 
kinematics of complex flows, the method of topological study of limiting 
streamlines is the most preferred one because the flow structure is best 
described by the topological properties inferred from computed streamlines. 
According to Legendre (1956), the study of the topology of three 
dimensional streamlines projected onto a no-slip body surface is 
computationally analogous to the experimental surface flow visualization. 
The skin-friction line field, according to the definition of Lighthill (1963), 
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can be regarded as the projection of the three-dimensional streamline field 
onto the body surface. This implies that the observed topological patterns of 
the limiting streamlines have significant implications for the flow evolution 
in the flow passage. 


It becomes quite convenient to analyze the kinematic aspects of computable 
limiting streamlines in terms of singular points of different kinds. So it 
becomes important to determine these singular points from which the flow 
structure can be sketched. In a set of velocity vectors, by definition the 
singular points are those, which have zero magnitude, but their directions are 
indeterminate. Around a critical point, the vector field varies according to 
the eigenvectors and the eigenvalues of the Jacobian of the field at the 
detected critical point. Depending on the divergence of the shear-stress 
vector and the Jacobian, singular points can be classified into two main 
topological types. These are saddle points and nodal points. The difference 
between these points can be observed very clearly by seeing the 
corresponding portraits (Figure 5.1). 

• Saddle points: - Through a saddle point only two sets of shear-stress lines 
pass. On each of these lines, the direction of the skin-friction vector field 

dv 

u — , u — 

_ oz oz 

saddle points, separation or attachment, both skin-friction directions are 
towards the critical point on one and both skin friction directions are 
away from this point on the other. In a two-dimensional flow there exists 


changes sign on a surface. Depending on the type of 
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Schematic diagrams of nodal and saddle points 

Figure 5.1 



special streamlines, called separation lines, which conveniently divide the 
flow into two distinct regions. Such streamlines originate at saddle points. 
The counterparts to these separation lines in three-dimensional flows are 
special streamsurfaces. This implies that the saddle points act as barriers in 
the field of vectors. 

• Nodal points:- These are points through which an infinite number of 
shear-stress lines pass. Nodal points are of two types, namely, regular 
nodal points (or nodes) and foci (or spiral nodes). Nodal points of 
attachment act as sources of skin-fnction lines that emanate from the 
point and spread out over the surface. Nodal points of separation, on the 
other hand, act as sinks where the skin friction lines circumscribing the 
body vanish. Foci and nodes differ in the orientation of their chosen field. 
For a singular point defined as a spiral node, all shear-stress lines spiral 
onto or out of this point. Emanating from this type of node are two 
straight critical lines. 

The motivation behind undertaking topological study of the vector field 
is to classify some mechanisms that may add to the heat-transfer 
capability of the system. A better understanding of the flow is essential. 
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5.2.1 Flow Topology on the Horizontal and the 
Vertical Planes of the Channel 


Figure 5.2 shows the streamline plot of the time-averaged flow on the 
horizontal midplane of the channel at z=H/2. The figure clearly shows the 
two symmetrical standing vortices. The flow separation and the wake 
formation are clearly visible. The recirculation in the wake region of the tube 
is clearly established. Figure 5.3 illustrates the streaklines, showing the path 
of the fluid particles on their way through the channel commencing at 
various starting points on a vertical plane. The flow pattern clearly shows the 
existance of a screw-like helical vortex motion in the wake region of the 
tube. Figure 5.4 shows the contour plot for the time averaged axial velocity, 

u, on the horizontal midplane of the channel at z=H/2. The recirculation 
region is located approximately at 0.39 diameters downstream away from 
the tube. The zero mean velocity line indicates the boundary of the 
recirculation zone. Figure 5.5 shows the contour plot of the transverse mean- 

velocity, V on the horizontal midplane of the channel. The profile is 
perfectly anti-Symmetric. The maximum value of transverse mean-velocity 
in the wake occurs at about 0.75 diameters downstream away from the tube. 

Figure 5.6 shows the distribution of time averaged pressure distribution 
along the tube wall. The pressure coefficient 

Cp = ( P - Pq + - yo i /? ), with Po the front stagnation 

2 2 
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Contour plot of time averaged axial velocity 
Figure 5.4 



Contour plot of mean transverse velocity 
Figure 5.5 





point pressure. The incomplete pressure recovery due to the recirculation 
behind the cylinder is clearly observed in the figure. 

Figure 5.7 shows the projection streamlines of the time averaged flow on 
the vertical midplane of the channel at y = B/2 . The wake of the cylinder 

exhibits a strong three-dimensional behavior. The strong normal velocity 
component, w is caused by the pressure gradient in the vertical direction in 
the wake region. This figure clearly shows the foot prints of the horseshoe 
vortices (marked as Si and Siin figure) and the reattachment line or the wake 
stagnation line (marked as S 3 in figure ). 

Figure 5.8 shows the projection of streamlines of the time-averaged flow on 
the transverse plane located at a distance of 1.2 time the radius of the tube 
from the center of the tube. A complex structure is discerned on the 
transverse plane behind the tube. Flow near the leading edge of the tube 
describes the deviation firom the two dimensional crossflow structure. 
Driven by the pressure gradient along the forward stagnation line of the tube, 
horseshoe vortices are formed in front of the tube. Two vortex systems, one 
at the junction of the top wall and tube and the other at the junction of the 
bottom wall and the tube, bend around the tube and form two sets of 
longitudinal vortex pairs. Kaul et al. (1985) describes the existence of the 
spiraling vortex systems, termed as the secondary vortex system, at the 
locations where the fluid is trapped between the end plates, tubes and the 
horseshoe vortices (also see Goldstein and Kami, 1984). As this secondary 
vortical flow filament on the bottom wall curls around the tube, it proceeds 
upward. The filament also begins to feel the effect of its counterpart coming 
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Streamlines on the vertical midplane of the channel 

Figure 5.7 





down from the top end wall. During their travel toward the symmetric 
midplane, the secondary vortex filaments interact with the horseshoe vortex 
systems and the vortex pair on both the walls get slightly shifted toward the 
symmetric midplane. 


5.2.2 Limiting Streamlines on the Tube Surface 
and the Bottom Plate 

The structure of the wall streamlines or limiting streamlines on the surface 
of the tube has investigated. This consists of the tangential components of 
velocity near the surface of the tube. The experimentalists have 
accomplished many investigations on the structure of limiting streamlines or 
the shear stresses lines on solid surfaces, but numerical computations are 
very few. Figure 5.9 shows such streamlines for the time-averaged flow. The 
mean flow field is symmetric so are the surface streamlines too. A 
bifurcation of the surface streamlines could be established. The bifurcation 
could be divided into two parts, a positive and negative bifurcation (Figure 
5.10 (a) and (b)). The positive bifurcation takes place along the forward 
stagnation line at the surface of the tube. This line is characterized by a 
single line on the incoming stream surface, which bifurcates into two as the 
positive bifurcation line (Figure 5.10 (a)). This phenomenon of positive 
bifurcation can also be seen in Figure 5.9. The left, middle and the right 
parts of the Figure 5.9 represents the forward stagnation line of the tube at 
0 = 0 and the line of reattachment at 6 = 7t . Figure 5.9, also reveals two 

other lines at an angle of OUn <0 < 0.8^ and 


77 



4.69 



Limiting streamlines on the tube surface 
Figure 5.9 




Positive and negative bifurcation lines 
Figure 5.10 


1.18;r < ^ < 1.3;zr respectively. These are the lines where the components of 
the shear stress vector in the direction of 0 are equal to 0. Homung and 
Perry (1984) call this a negative stream surface bifurcation, which can also 
be seen in Figure 5.10 (b). Two streamlines combine to form a single 
streamline along the negative bifurcation line. These lines can be termed as 
separation lines along which the boundary layer separates from the .tube 
surface. Figure 5.9 makes it obvious that the separation angle is small in the 
vicinity of the midplane of the tube and becomes larger in the vicinity of end 
walls, because of the greater flow speed in the plane of sy mm etry of the 
tube. 

Figure 5.11 gives the streamlines of the time-averaged flow in the region 
close to the plate. A great deal of information about the flow and heat 
transfer characteristics can be extracted from this figure. A saddle point of 
separation and a horseshoe vortex system is observed from the figure. In 
Figure 5.11, the incoming flow does not separate in the traditional sense but 
reaches a stagnation or saddle point of separation (marked by A on the 
figure) and goes around the body. The nodal point of attachment (marked as 
C) and the separation lines which form circular arcs across the tube are 
clearly visible in the figure. The flow above the lower wall hits the front of 
the tube. A significant part of it moves downward and creates a region of 
reversed flow in front of the stagnation line. On each side of the tube, one 
finds a region of converging streamlines (marked as G). These are the traces 
of horseshoe vortex. Behind the body one finds two areas of swirling flow 
(marked as E) which are the footprints of an arch vortex. Finally, there is a 
wake stagnation point (marked F) further downstream of the body. 
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5.3 Heat Transfer Characteristics 


In this section a detailed study of heat transfer characteristics for the fin tube 
element (Figure 1.3) is presented. An effort has been made to understand the 
heat transfer characteristics on the basis of the flow pattern for the fm-tube 
heat exchanger module, which has been discussed in detail in previous 
section. 

5.3.1 Temperature Distribution 

Figure 5.12 shows the time-averaged temperature contours close to the 
bottom plate. Figure 5.13 shows the corresponding surface plot. The 
temperature field shows a dip near the forward stagnation line. This is due to 
the entertainment of cold fluid from the horizontal midplane of the channel. 
A high temperature in the wake region is also observed. In the wake region, 
the fluid recirculates with low velocity. Hence it gets heated more as 
compared to the neighboring regions. A low temperature is observed in a 
region wrapping the circumference of the tube. This is due to creation of the 
horseshoe vortices, which constitute the high-velocity spiraling motion of 
the fluid that sweeps cooler fluid from the outer regions of the thermal 
boundary layer toward the wall and consequently a low temperature is 
observed in this region. 

Figure 5.14 shows the time-averaged temperature contours on the vertical 
midplane of the channel. The figure clearly reveals the thermal boundary 
layer development along the bottom plate. The contours, as expected, are 


82 



Time-averaged temperature contours close to the bottom plate 

Figure 5.12 
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ontours on the bottom plate 





Time-averaged temperature contours on the vertical midplane of the 

channel 
Figure 5.14 
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Temperature contours on the vertical midplane 


symmetric with respect to z=]-I/2. The gradients arc clustered primarily in 
regions in front of the tube. 


5.3.2 Nusselt Number and the Performance of the 
Heat Exchanger Module 


Figure 5.15 shows the contour map of local Nusselt number obtained from 
the time-averaged temperature field for the bottom plate. Several interesting 
features are observed from the figure. Primarily foil wing three regions may 
be distinguished. At the leading edge of the bottom plate, Nusselt number 
reaches its maximum and decreases gradually. At the leading edge, the 
cooler fluid comes in contact with the hot solid wall for the first time; hence 
the heat transfer is maximum. The gradual decrease in the Nusselt number is 
attributed due to the boundary layer development on the channel wall. This 
observation does not match with the experimental results of O’Brien et al. 
(2000) for the same Reynolds number. This is because in the experiment a 
flow development region was provided in the upstream of the test section. 
Consequently, in the experiment, the velocity boundary layer and the 
thermal boundaiy layer is approximately fully developed as it enters the test 
section. Present numerical simulation does not consider any flow 
development section. The uniform velocity and temperature profiles have 
been deployed in this study. 

An abmpt increase of Nusselt number in front of the tube is observed. This 

results from the formation of horseshoe vortex system that consists of tw'o 
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Contour map of local Nusselt number close to bottom plate 

Figure 5.15 





counter-rotating longitudinal vortices. As the fluid approaches the 
stagnation line of the circular tube, it slows down and its pressure increases. 
The smaller velocity in the boundary layer in the vicinity of the bottom 
plate, which supports the circular tube, leads to a smaller pressure increase. 
Thus, the induced pressure gradient causes the flow towards the bottom 
wall that interacts with the main stream. The fluid rolls up forming vortices, 
which finally wraps around the front half of the tube and extends to the rear 
of the tube, as shown schematically in Figure 5.16 (after Goldstein and 
Kami (1984). The spiraling motion of the horseshoe vortices brings about a 
better mixing and the heat transfer in this region is enhanced significantly. 
This is in good qualitative agreement with the experimental results presented 
by O’Brien et al. (2000). The quantitative comparisons have not been made 
because in experiments, the tube and the wall temperatures are not constant 
since both convective and conductive heat transfer take place from the tube 
and the wall. Present numerical simulation assumes the tube wall and fins to 
be at constant temperature. In addition, in the experiments, the sidewalls are 
no-slip walls whereas our simulation considers the sidewalls as the free-slip 
walls. 

The Nusselt number is low in the wake region as shown in Figure 5.15. The 
poor heat transfer in this region is attributed to the separated dead water zone 
with fluid recirculating at a low velocity. This trend also matches 
qualitatively quite well with the experimental result of O’ Brien et al. 
( 2000 ). 
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Three-dimensional boundary layer separation and horseshoe vortex 
system in the region of interaction between mainstream boundary 
layer and tube protruding from the wall( Goldstien and Kami (1984)) 

Figure 5.16 


Figure 5.17 shows the distribution of span-averaged Nusselt number in the 
streamwise direction. At the leading edge of the fin, Nusselt number reaches 
its maximum and decreases gradually. There is an abrupt increase in 
Nusselt number in front of the tube. The reasons for this behavior are 
explained in this section. 

Figure 5.18 represents the transverse variation in the local Nusselt number 
in the wake region at a distance of x/R = 1.2 from the center of the tube. 
The two peaks seen in the figure result from the downstream longitudinal 
bportion of the horseshoe vortices. The low heat transfer in the wake region 
(-0.5< y/R <0.5) is clearly established from the figure. 
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Span-averaged Nusselt number distribution in the streamwise 

direction 
Figure 5.17 



Chapter 6 

Conclusion and Scope for Future 
Work 

6.1 Conclusions 

A three-dimensional numerical study on the flow and heat transfer 
characteristics in a narrow rectangular duct with a built-in circular tube in 
cross-flow has been performed. The duct was designed to simulate a 
passage, formed by two neighboring fins in a fin-tube heat exchanger. The 
governing conservation equations of fluid flow are derived in its integral 
form. An explicit control volume formulation devised by Eswaran and 
Prakash (1998) is used to solve the equations. 

The flow field is completely different from the two dimensional flows. The 
third velocity component w is surprisingly very large in the front and the 
near wake of the tube. In the rear of the tube it leads in accompaniment with 
the back flow to a screw-like motion of a helical vortex tube. Limiting 
streamlines on the tube surface have shown the typical separation lines. 
Limiting streamlines on the bottom wall have clearly established the 
presence saddle point of separation, the nodal point of attachment in front of 
the tube-bottom wall junction and the horseshoe vortex system that wraps 
around the tube and extend to the rear of the tube. The separation line 
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emanating from the saddle point is also clearly established. The local span- 
averaged Nusselt number and the iso-Nusselt number distribution clearly 
establish the high heat transfer near the leading edge of the fins and in the 
region influenced by the horseshoe vortex. The poor heat transfer in the dead 
water zone is observed. A double peak in local fm-surface heat transfer in 
the stagnation region is clearly observed. Although no quantitative 
comparisons can be made with the available experimental results of O’Brien 
(2000) because of the different flow and geometric conditions, the 
qualitative comparisons are meaningful and it matches quite well with 
experimental results. 

6.2 Scope for Future Work 

The results of this work validate three-dimensional unsteady flow models as 
a useful technique in improving the heat transfer efficiency for air-cooled 
heat exchangers. As the results reveal, regions of relatively high heat transfer 
are constrained to small regions just upstream and to the sides of the tube 
and in the immediate vicinity of the horseshoe vortices downstream of the 
tube. The heat transfer performance can be improved by placing a delta 
winglet pair behind the tube. The winglet pair can be mounted following two 
different arrangements. The conventional way of placing the winglet pair is 
such that the transverse distance between the leading edges is less than the 
transverse distance between the trailing edges. Since the vortical structure 
due to such an arrangement shows common flow in the direction of the wall 
(bottom plate of the channel), the configuration is called ‘common flow 
down’. However, Torii et al. (2000) have suggested a special arrangement ol 


94 



placing the winglets behind the circular tube. In this arrangement, the 
distance between the trailing edges of the winglet pair is less than the 
distance between the leading edges. The vortical structure entails common 
flow away from the wall in such an arrangement. Hence the configuration is 
called common flow up’. The salient feature of this arrangement is creation 
of a virtual constriction of the flow passage near the rear stagnation point of 
the tube. This constriction brings about considerable separation delay and 
reduces form drag. Due to the accelerating stream, the poor heat transfer 
zone in the wake of the cylinder is removed. The future investigations with 
the winglet pair should include both the configurations, such as ‘common 
flow down’ and common flow up’. Biswas and Mitra (1988) have 
conducted comparative study using delta winglet pair and the rectangular 
winglet pair as the vortex generators. As a future computational study both 
delta winglet pair and rectangular winglet pair can be used and the condition 
for the optimum heat transfer can be determined. 
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Appendix 


Attachment and separation in three-dimensional flows 


At a very small normal distance z from the solid surface, 

^du u dv 

so 

V oz z oz z 


II 


(A-1) 

^ (du 

“4 

T 


where 

^ uZ C/Z j 

^ 

M 

(A-2) 

where, represents the skin friction vector field. 


From (A-1 ) we have. 



u = e^z, 

for solenoidal flow field 

div^ = 0 where q = (m,v,w) 

II 

((A-3) 


or, — = —A z , where A = div 
oz 

On integrating we have 

W = — Az" (A-4) 

2 

from equation (A-4) it is obvious that streamlines very close to the 

-> 

surface will be parallel to the surface (w=0) provided 

Critical Points: Critical points are those where 6:^ = 0 or the velocity 
vector vanishes, therefore its directions are indeterminate. 

Separation Points: are those points where the vertical (normal) 
component of velocity w>0 



If w > Ofrom equation (A-4) we have A < 0 , i.e . the divergence of skin 
friction vector is negative at the separation point. 

Attachment Point: are those points where the normal component of 
velocity w<0 

w<0, or A>0 (A-5) 

The divergence of skin friction vector is positive at the attachment point. 

Skin Friction Lines : are curves such that tangent at any point give 
direction of the or the skin friction vector field at that point. 

From equation (A-1) it is obvious that for the location very close to the 
surface, streamlines lie closely along with the skin friction lines. 

Critical Points can be divided primarily of two types, Saddle points and 
nodal points . The nature of the critical point depends on the Jacobian of 

the shear stress vector (^^)and the divergence of the skin friction 
vector field. 

Saddle Points: For any initial point in the neighborhood of this critical 
point, with passage of time, there will be attraction towards this point 
along one direction and repulsion from this point along another 
direction. 


Mathematically 


r 


J 




y) 




J (>% »yo ) 



<0 


a»(j£o>yo) 


or it corresponds to one positive or negative eigen values. If A < Oand 
J < 0 => saddle point of separation 

if A > 0 and J < 0 saddle point of attachment 



Saddle 


Saddle 


Saddle point of separation Saddle point of reattachment 

Nodal Points: are the points through which an infinite number of shear 
stress lines pass. Nodal Points arp of two types, regular nodal points (or 
nodes) and foci (or spiral nodes). 

Nodes: Mathematically J > 0 for the nodal points 




J >0,A>0,w<0 J >0,A<0,>v>0 

The above two are regular nodes. For the stable node (sink), both the 
eigen values are negative. For the unstable node, both the eigen values 
are positive. 

Spiral node: Here the skin friction lines spiral onto or out of this point. 
This correspond to complex eigen values. If the real part of the eigen 
value Is negative It is a stable spiral and if the real part of the eigen 
value is positive, it corresponds to an unstable spiral. 




Stable Spiral 


Unstable Spiral 



